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Abstract 



We present results of lattice QCD simulations with mass-degenerate up and down 
and mass-split strange and charm (N{ = 2 + 1 + 1) dynamical quarks using Wilson 
twisted mass fermions at maximal twist. The tuning of the strange and charm 
quark masses is performed at two values of the lattice spacing a ~ 0.078 fm and 
a ~ 0.086 fm with lattice sizes ranging from L ~ 1.9 fm to L ~ 2.8 fm. We measure 
with high statistical precision the light pseudoscalar mass mps and decay constant 
fps in a range 270 < tops < 510 MeV and determine the low energy parameters /o 
and Z3 4 of SU(2) chiral perturbation theory. We use the two values of the lattice 
spacing, several lattice sizes as well as different values of the light, strange and 
charm quark masses to explore the systematic effects. A first study of discretisation 
effects in light-quark observables and a comparison to N{ = 2 results are performed. 
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1 Introduction and Main Results 



The beginning of this century has assisted to radical improvements in theory, 
algorithms and supercomputer technology, leading to a far increased ability to 
solve non-perturbative aspects of gauge field theories in a lattice regularised 
framework. Following this path of improving the lattice setup, in this paper, we 
are reporting about our experiences and results when considering in addition 
to the u, d light dynamical flavours also the effects of the strange and charm 
sea quarks. By including a dynamical charm, we are now able to directly study 
its contribution to physical observables and to quantify the so far uncontrolled 
systematic effect present in lattice QCD simulations where the charm flavour 
in the sea is absent. 

A number of different lattice fermion formulations are being used by several 
lattice groups, see refs. [1, 2] for recent reviews. Here, we adopt a particular 
type of Wilson fermions, known as the Wilson twisted mass formulation of 
lattice QCD (tmLQCD), introduced in [3,4]. This approach is by now well 
established, with many physical results obtained with two light degenerate 
twisted mass flavours (JVf = 2) by our European Twisted Mass (ETM) Col- 
laboration, see refs. [5-22]. For a review see ref. [23]. In the tmLQCD formula- 
tion a twisted mass term is added to the standard, unimproved Wilson-Dirac 
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operator, and the formulation becomes especially interesting when the theory 
is tuned to maximal twist [4]. The major advantage of the lattice theory tuned 
to maximal twist is the automatic 0(a) improvement of physical observables, 
independently of the specific type of operator considered, implying that no 
additional, operator specific improvement coefficients need to be computed. 
Other advantages worth to mention are that the twisted mass term acts as an 
infrared regulator of the theory and that mixing patterns in the renormalisa- 
tion procedure are expected to be simplified. 

Detailed studies of the continuum-limit scaling in the quenched approxima- 
tion [24-27] and with two dynamical quarks [7,10,17,28] have demonstrated 
that, after an appropriate tuning procedure to maximal twist, lattice artefacts 
not only follow the expected 0(a 2 ) scaling behaviour [4], but also that the re- 
maining 0(a 2 ) effects are small, in agreement with the conclusions drawn in 
ref. [29]. 

The only exception seen so far is the neutral pseudoscalar mass, which shows 
significant 0(a 2 ) effects. This arises from the explicit breaking of both parity 
and isospin symmetry, which are however restored in the continuum limit 
with a rate of 0(a 2 ) as shown in [4] and numerically confirmed in refs. [17,30]. 
Moreover, a recent analysis suggests that isospin breaking effects strongly 
affect only a limited set of observables, namely the neutral pion mass and 
kinematically related quantities [31,32]. 

In this paper we report on simulations with twisted mass dynamical up, down, 
strange and charm quarks. We realise this by adding a heavy mass-split dou- 
blet (c, s) to the light degenerate mass doublet {u,d), referring to this setup 
as Nf = 2 + 1 + 1 simulations. This formulation was introduced in [33,34] and 
first explored in [35]. As for the mass-degenerate case, the use of lattice action 
symmetries allows to prove the automatic 0(a) improvement of physical ob- 
servables in the non-degenerate case [33, 34] . First accounts of our work were 
presented at recent conferences [36,37]. Recently, results with Nf — 2 + 1 + 1 
staggered fermions have been reported in [38-40]. The inclusion of the strange 
and charm degrees of freedom allows for a most complete description of light 
hadron physics and eventually opens the way to explore effects of a dynamical 
charm in genuinely strong interaction processes and in weak matrix elements. 

Here, we concentrate on results in the light-quark sector using the charged 
pseudoscalar mass m PS and decay constant fps as basic observables involving 
up and down valence quarks only. In fig. 1 we show the dependence of (a) 
mp S /2Bofii and (b) fps as a function of the mass parameter 2B fii, together 
with a fit to SU(2) chiral perturbation theory (xPT) at the smallest value of 
the lattice spacing of a ~ 0.078 fm and lattice gauge coupling (3 = 1.95. We 
summarise the fit results for the low energy constants in table 1. These are 
the main results of this paper. 
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(a) (b) 

Fig. 1. (a) The charged pseudoscalar mass ratio mf> s /(2.E?o/i;) an d (b) the pseu- 
doscalar decay constant /ps as a function of 2B§[i\ fitted to SU(2) chiral perturba- 
tion theory, see table 1. The scale is set by the value of 2Bqhi at which the ratio 
fps~°°^ I ' m PS~°°^ assumes its physical value [41] f^/m^ = 130.4(2)/135.0 (black star). 
The lattice gauge coupling is f3 = 1.95 and the twisted light quark mass ranges from 
am = 0.0025 to 0.0085, see eq. (3) for its definition, corresponding to a range of 
the pseudoscalar mass 270 < ?nps < 490 MeV. The kaon and D meson masses are 
tuned to their physical value, see table 4. The lightest point (open symbol) has not 
been included in the chiral fit, see the discussion in section 3.2. 



A comparison between data obtained with Nf = 2 + 1 + 1 and Nf = 2 flavours 
of quarks - see sections 3.4 and 4, and ref. [17] - reveals a remarkable agreement 
for the results involving light-quark observables such as the pseudoscalar mass 
and decay constant or the nucleon mass. This provides a strong indication in 
favour of the good quality of our data in this new setup. In particular, barring 
cancellations due to lattice discretisation errors, these results would suggest 
that the dynamical strange and charm degrees of freedom do not induce large 
effects in these light-quark observables. In the Nf = 2 case, data collected 
at four values of the lattice spacing have allowed us to properly quantify all 
systematic errors present in the determination of light-quark observables [17]. 
In this first work with Nf = 2 + 1 + 1 flavours, we consider data at two close 
values of the lattice spacing, while we defer to a forthcoming publication the 
inclusion of additional ensembles at a significantly lower lattice spacing and a 
more complete analysis of the systematic effects. 

The rest of this paper is organised as follows. In section 2 we describe the 
gauge action and the twisted mass fermionic action for the light and heavy 
sectors of the theory. The realisation of 0(a) improvement at maximal twist 
is also presented. In section 3 we define the simulation parameters, describe 
the tuning to maximal twist as well as the tuning of the strange and charm 
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P = 1.95 

h 3.70(7X26)" 
k 4.67(3)(10) 
fo [MeV] 121.14(8)(19) 
U/fo 1.076(2) (2) 

2Bofi u ,d/K 1.032(21)(3) 

(r 2 )^ LO [fm 2 ] 0.724(5)(23) 

r$/a(P = 1.95) 5.71(4) 
r*(/3 = 1.95) [fm] 0.447(5) 
gQ3 = 1.95) [fm] 0.0782(6) 

Table 1 

Results of the fits to SU(2) xPT for the ensemble at /3 = 1.95. Predicted quantities 
are: the low energy constants £3,4, the charged pseudoscalar decay constant in the 
chiral limit fo, the mass ratio 2Bo/J,i/mp S at the physical point and the pion scalar 
radius (r 2 )^ L °. The first quoted error is from the chiral fit at (3 = 1.95, the second 
error is the systematic uncertainty that conservatively accommodates the best fitted 
central values of the three fits reported in table 9, section 4. The small error on the 
quoted lattice spacing comes exclusively from the fit at /3 = 1.95. The scale is 

set by fixing the ratio /p^ _ °°' / m ps~°°' = fir/ m -K = 130.4(2)/135.0 to its physical 
value [41]. The chirally extrapolated Sommer scale is determined separately and 
not included in the xPT fits. For a comparison with the N{ = 2 ETMC results, 
see [17]. 



quark masses and the relevance of discretisation effects. Section 4 includes a 
discussion of the fits to SU(2) %PT also for data on a slightly coarser lattice, 
a PS 0.086 fm, and provides a first account of systematic uncertainties. Our 
conclusions and future prospects are summarised in section 5. 



2 Lattice Action 

The complete lattice action can be written as 

S = S g + S l + S h , (1) 



where S g is the pure gauge action, in our case the so-called Iwasaki action [42, 
43], S[ is the twisted mass Wilson action for the light doublet [3,4] and Sh the 
one for the heavy doublet [33,34]. 
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2.1 Gauge action 



The Iwasaki gauge action [42,43] includes besides the plaquette term 
also rectangular (1 x 2) Wilson loops 

S 9 = f E f &0 E {1 - Re TriU^+h E {1 - Re Tr([/^)} ) , (2) 

with /3 = 6/(7q the bare inverse coupling, &! = —0.331 and the normalisation 
condition b = 1 — 8&i. 

The choice of the gauge action is motivated by the non trivial phase struc- 
ture of Wilson-type fermions at finite values of the lattice spacing. The phase 
structure of the theory has been extensively studied analytically, by means 
of chiral perturbation theory [44-50], and numerically [51-56]. These studies 
provided evidence for a first order phase transition close to the chiral point for 
coarse lattices. This implies that simulations at non-vanishing lattice spacing 
cannot be performed with pseudoscalar masses below a minimal critical value. 

The strength of the phase transition has been found [53, 56] to be highly 
sensitive to the value of the parameter b\ in the gauge action in eq. (2). 
Moreover, in [35] it was observed that its strength grows when increasing 
the number of flavours in the sea from Nf — 2 to Nf — 2 + 1 + 1, at otherwise 
fixed physical situation. Numerical studies with our Nf = 2 + 1 + 1 setup have 
shown that the Iwasaki gauge action, with b\ = —0.331, provides a smoother 
dependence of phase transition sensitive quantities on the bare quark mass 
than the tree- level- improved Symanzik [57,58] gauge action, with b\ = —1/12, 
chosen for our Nf = 2 simulations. 

Another way to weaken the strength of the phase transition is to modify the 
covariant derivative in the fermion action by smearing the gauge fields. While 
the main results of this work do not use smearing of the gauge fields, we 
report in section 3.7 on our experience when applying a stout smearing [59] 
procedure, see also [60]. 

2.2 Action for the Light Doublet 

The lattice action for the mass degenerate light doublet (u, d) in the so called 
twisted basis reads [3,4] 

s i = « 4 Efe(X) \ D [U] +m ,i + im^}Xi(x)} , (3) 

X 
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where mo,i is the untwisted bare quark mass, [i\ is the bare twisted light quark 
mass, r 3 is the third Pauli matrix acting in flavour space and 



Dp] = \ [ 7 „ (v„ + v;) - av;v M 



2 

is the massless Wilson-Dirac operator. V M and V* are the forward and back- 
ward gauge covariant difference operators, respectively. Twisted mass light 
fermions are said to be at maximal twist if the bare untwisted mass m ,z is 
tuned to its critical value, m crit , the situation we shall reproduce in our sim- 
ulations. The quark doublet xi — (Xu, Xd) i n the twisted basis is related by a 
chiral rotation to the quark doublet in the physical basis 

^fys = e i^ Xu ^Hys = ^n B r 3 ? (4) 

where the twisting angle oji takes the value \u>i\ — > | as |m 0j z — m cr i t | — > 0. We 
shall use the twisted basis throughout this paper. 



2.3 Action for the Heavy Doublet 



We introduce a dynamical strange quark by adding a twisted heavy mass- 
split doublet Xh — (XcXs), thus also introducing a dynamical charm in our 
framework. As shown in [34], a real quark determinant can in this case be 
obtained if the mass splitting is taken to be orthogonal in isospin space to the 
twist direction. We thus choose the construction [33,34] 

s h = a 4 J2 iXh(x) [D[U] + m 0;h + ifi^n + /i«5T 3 ] Xh(x)} , (5) 

X 

where m ,/i is the untwisted bare quark mass for the heavy doublet, jj, a the 
bare twisted mass - the twist is this time along the T\ direction - and /ig the 
mass splitting along the T3 direction. 

The bare mass parameters \x a and /i«s of the non-degenerate heavy doublet are 
related to the physical renormalised strange and charm quark masses via [33] 

(m s ) R = Zp 1 (fj, a - Zp/Z s ps) , /gx 
(m c ) R = Zp 1 (fj, a + Zp/Z s us), 

where Zp and Z$ are the renormalisation constants of the pseudoscalar and 
scalar quark densities, respectively, computed in the massless standard Wilson 
theory. 

A chiral rotation analogous to the one in the light sector transforms the heavy 
quark doublet from the twisted to the physical basis 
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where the twisting angle Uh takes the value \u>h\ — > § as |m 0j /i — m crit | — > 0. 
2.4 0(a) improvement at maximal twist 



One of the main advantages of Wilson twisted mass fermions is that by tuning 
the untwisted bare quark mass to its critical value, automatic 0(a) improve- 
ment of physical observables can be achieved. 

Tuning the complete Nf = 2 + 1 + 1 action to maximal twist can in principle be 
performed by independently choosing the bare masses of the light and heavy 
sectors amoj and am ,/i, resulting, however, in a quite demanding procedure. 
On the other hand, properties of the Wilson twisted mass formulation allow 
for a rather economical, while accurate alternative [4,34,35], where the choice 
am 0t i = am 0j h = 1/2/c — 4 is made, and the hopping parameter k has been 
introduced. 

Tuning to maximal twist, i.e. n = K crit , is then achieved by choosing a parity 
odd operator O and determine am cr i t (equivalently K C ru) such that O has van- 
ishing expectation value. One appropriate quantity is the PCAC light quark 
mass [29,52,53] 

E x (^(x,t)P*(0)) 

-pcac= 2Ex(if(x>t)if(0)> , « = 1,2, (8) 

where 

A^ix) = xi(xh^xi(x) , P t a (x) = jato-y^xiix) , (9) 

and we demand mpcAc = 0. For the quenched [25] and the Nf — 2 case [17], 
this method has been found to be successful in providing the expected 0(a) 
improvement and effectively reducing residual O(o?) discretisation effects in 
the region of small quark masses [29] . 

The numerical precision required for the tuning of m PC Ac to zero has been 
discussed in [8]. Contrary to the Nf — 2 case [5,8], where this tuning was 
performed once at the minimal value of the twisted light mass considered in 
the simulations, we now perform the tuning at each value of the twisted light 
quark mass /ii and the heavy-doublet quark mass parameters \i a and fis- This 
obviously leaves more freedom in the choice of light quark masses for future 
computations. 

Although theoretical arguments tell us that 0(a) improvement is at work in 
our setup, a dedicated continuum scaling study is always required to accurately 
quantify the actual magnitude of 0(a 2 ) effects. In section 3.4 we provide a first 
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indication that such effects are indeed small, at least for the here considered 
light meson sector; currently ongoing computations at a significantly smaller 
lattice spacing will allow for a continuum limit scaling analysis in this setup. 



3 Simulation Details 

3.1 Simulation Ensembles 

We performed simulations at two values of the lattice gauge coupling /3 = 1.90 
and 1.95, corresponding to values of the lattice spacing a 0.086 fm and 
a rH 0.078 fm, respectively. The parameters of each ensemble are reported 
in table 2. The charged pion mass mpg ranges from 270 MeV to 510 MeV. 
Simulated volumes correspond to values of mpgL ranging from 3.0 to 5.8, 
where the smaller volumes served to estimate finite volume effects, see table 3. 
Physical spatial volumes range from (1.9 fm) 3 to (2.8 fm) 3 . 

As already mentioned, the tuning to K, cr u was performed independently for 
each value of the mass parameters afii, afi a and afi$. The mass parameters of 
the heavy doublet a/v and afis reported in table 2 are related to the strange 
and charm quark masses. In particular, they are fixed by requiring the simu- 
lated kaon and D meson masses to approximately take their physical values, 
as discussed in section 3.3. The simulation algorithm used to generate the 
ensembles includes in the light sector, a Hybrid Monte Carlo algorithm with 
multiple time scales and mass preconditioning, described in ref. [61], while 
in the strange-charm sector a polynomial hybrid Monte Carlo (PHMC) algo- 
rithm [62-64]; the implementation of ref. [65] is publicly available. 

The positivity of the determinant of the Dirac operator is a property of the 
mass-degenerate Wilson twisted mass action, which does not necessarily hold 
in the non degenerate case for generic values of the mass parameters \x a and 
^.[MThe positivity is monitored by measuring the smallest eigenvalue Ah, m in 
of Q^Qhi where Qh = l5T%D h and D h is the Wilson Dirac operator of the non- 
degenerate twisted mass action in eq. (5). We observe that Ah, m in is roughly 
proportional to the renormalised strange quark mass squared. Since we choose 
the mass parameters \x a and fis such that the strange quark takes its physical 
value, a spectral gap in the distribution of Q^Qh is observed, implying that 
the determinant of does not change sign during the simulation. While this 
is sufficient for the purpose of this study, we shall provide a detailed discussion 
of this issue in a forthcoming publication. 

1 Notice however that the positivity of the determinant is guaranteed for ^ > 
A [33,34]. 
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Ensemble 





^crit 




a/j, a 




(L/a) 3 x T/a 


A ^9 

AOU.OZ 


i on 
i .yu 


U. IDoZ / zu 


u.uuou 


U. 10U 


n i an 
u. iyu 


"39 3 v P.A 


A /in ^9 




U. / uu 


n nn/in 

U.UU4U 






^9 3 v R/l 


A 4D 94 




U. IDoZ / uu 


U.UU4LU 






94 3 y /IQ 


A40.20 




n 1632700 


0.0040 






20 3 x 48 


A 50 39 




1639670 


0050 






39 3 x 64 


A fift 94 

nUU. Zi4 






OOfiO 






94 3 v 4R 


A 80 94 






0080 






94 3 v 48 


a i nn 94 




U. lUOiUuU 


01 00 






94 3 v 48 


A1UU.Z4S 




u. looiyou 


n ni nn 
U.U1UU 




u.iy / 


Z4 X 4o 


B25.32 


1.95 


0.1612420 


0.0025 


0.135 


0.170 


32 3 x 64 


B35.32 




0.1612400 


0.0035 






32 3 x 64 


B55.32 




0.1612360 


0.0055 






32 3 x 64 


B75.32 




0.1612320 


0.0075 






32 3 x 64 


B85.24 




0.1612312 


0.0085 






24 3 x 48 



Table 2 

Summary of the Nf = 2 + 1 + 1 ensembles generated by ETMC at two values of the 
lattice coupling = 1.90 and = 1.95. From left to right, we quote the ensemble 
name, the value of inverse coupling 0, the estimate of the critical value K cr it, the light 
twisted mass a/ii, the heavy doublet mass parameters a/j, a and a/j,g and the volume 
in units of the lattice spacing. Our notation for the ensemble names corresponds to 
X.fii.L, with X referring to the value of used. The run Al00.24s is used to control 
the tuning of the strange and charm quark masses. 



Ensemble 


rapCAc/w 


mp S L 


Tint((P» 


ri nt (amps) 


r in t(ampcAc) 


A30.32 


-0.123(87) 


3.97 


118(55) 


2.7(4) 


46(19) 


A40.32 


-0.055(55) 


4.53 


103(48) 


4.1(7) 


51(21) 


A40.24 


-0.148(83) 


3.48 


132(57) 


< 2 


35(12) 


A40.20 


-0.051(91) 


2.97 


55(25) 


2.9(7) 


26(12) 


A50.32 


0.064(24) 


5.05 


50(19) 


3.0(5) 


21(7) 


A60.24 


-0.037(50) 


4.15 


28(8) 


2.0(2) 


13(4) 


A80.24 


0.020(19) 


4.77 


23(7) 


2.4(3) 


10(2) 


A100.24 


0.025(18) 


5.35 


18(5) 


2.3(3) 


13(3) 


A100.24s 


0.045(18) 


5.31 


18(5) 


6.2(1.1) 


18(5) 


B25.32 


-0.185(69) 


3.42 


65(25) 


3.6(6) 


26(9) 


B35.32 


0.009(34) 


4.03 


54(19) 


5.5(8) 


41(14) 


B55.32 


-0.069(13) 


4.97 


12(3) 


< 2 


8(2) 


B75.32 


-0.047(12) 


5.77 


14(4) 


3.3(5) 


13(3) 


B85.24 


-0.001(16) 


4.66 


15(4) 


2.2(2) 


11(2) 



Table 3 

For each ensemble, from left to right the values of mpcAc/rt, m PsL, the integrated 
autocorrelation time of the plaquette, mps and mpcAC m units of the trajectory 
length. Every ensemble contains 5000 thermalised trajectories of length r = 1, 
except A40.24 which contains 8000 trajectories. 

To generate correlators we use stochastic sources and improve the signal-to- 
noise ratio by using the "one-end trick" , following the techniques also employed 
in our N { = 2 simulations [8] . We have constructed all meson correlators with 
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,3=1.90 • L/a = 24 
• 13= 1.95 ▲ L/a = 32 



I 



i 



0.15 0.20 

2B„ W [GcV 2 ] 



0.25 



Fig. 2. The ratio mpcAc/w for the ensembles at f3 = 1.90 and 1.95 at the largest sim- 
ulated volumes and as a function of 2Bq^i. For both ensembles the ratio mpcAC/Vz 
satisfies the 10% level criterion, except for the lightest point at f3 = 1.90 and 
f3 = 1.95 (open symbols), also affected by larger statistical errors. We assume 
Za = 1, while the actual value Za < 1 can only improve all tuning conditions. 

local (L), fuzzed (F) and Gaussian smeared (S) sources and sinks. The use of 
smeared or fuzzed sources has stronger impact on the extraction of the kaon 
and D meson masses; results for the latter are reported in section 3.3, while 
a companion paper [66] discusses the adopted strategy for the less straight- 
forward determination of these masses in the unitary iV f = 2 + 1 + 1 Wilson 
twisted mass formalism. 



3.2 Tuning to Maximal Twist 



To guarantee 0(a) improvement of all physical observables while also avoiding 
residual 0(a 2 ) effects with decreasing pion mass, the numerical precision of 
the tuning to maximal twist - quantified by the deviation from zero of mpcAc 
- has to satisfy |i^mpcAc/ /^U,M<r,w ~ o,A QC d [5,8,17]. The left-hand side 
contains the renormalised ratio of the untwisted mass over the twisted light- 
quark mass. A similar condition should be fulfilled by the error on this ratio. 
For the current lattice spacings, oKqcd ~ 0.1, while the values of the axial 
current renormalisation factor have not yet been determined. Nevertheless, 
since Z& enters as an 0(1) multiplicative prefactor, and it is expected to be 
Zj^ < 1 for our ensembles^} we adopt the conservative choice Za = 1 in 
verifying the tuning condition. 



2 Preliminary determinations of Za from ongoing dedicated runs with four degen- 
erate light flavours, indicate that Za ~ 0.7 — 0.8 for the ensembles considered in 
this work. 
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Satisfying this constraint clearly requires a good statistical accuracy in the 
determination of the PCAC mass. The values of 7Bpcac/ l^i reported in table 3 
and shown in fig. 2 are well satisfying the tuning condition to maximal twist, 
with the exception of the lightest mass point at /3 — 1.90 and (5 = 1.95. We 
notice that the autocorrelation time of mpcAc reported in table 3 grows with 
decreasing values of the light quark mass fii, thus rendering the tuning more 
costly for the two lightest points. For the ensemble B25.32, we are currently 
performing a new simulation aiming at a more accurate tuning to K CTlt . We 
are also testing a reweighting procedure [36] in k on the same ensemble, in 
view of applying it to the other not optimally tuned ensemble A30.32, and 
to future simulations. In what follows, we use the lightest mass points for 
consistency checks, and we exclude them from the final %PT fits. We also 
remind the reader that the small deviations from zero of ampcAc w iU only 
affect the 0(a 2 ) lattice discretisation errors of physical observables [8]. 



3. 3 Tuning of the Strange and Charm Quark Masses 

The mass parameters fi a and fi$ in the heavy doublet of the action in eq. (5) 
can in principle be adjusted so as to match the renormalised strange and charm 
quark masses by use of eq. (6). In practise, in this work, we fix the values of 
jjjfj and fis by requiring that the simulated kaon mass rriK and D meson mass 
mo approximately take their physical values. 

A detailed description of the determination of the kaon and D meson masses is 
separately given in [66], while figures 3(a) and 3(b) show the resulting depen- 
dence of (2m 2 K — m 2 PS ) and mp upon the light pseudoscalar mass squared for 
both ensembles, and compared with the physical point. Table 4 summarises 
their numerical values, while the corresponding values for a\i a and afis are 
given in table 2. Observe also that, in order to be able to properly tune the 
strange and charm quark masses to their physical values, a\i a must be chosen 
larger than a/is, since (see eq. (6)) the ratio Zp/Zs is significantly smaller 
than one [66]. 

While the kaon and D meson masses at (3 = 1.95 are sufficiently well tuned 
to their physical values, the ensembles at (5 = 1.90 with a^s = 0.190 carry a 
heavier kaon mass. The latter is instead visibly closer to its physical value for 
a^s = 0.197, as can be inferred from figure 3(a). We are currently performing 
simulations with a/is = 0.197 for other light quark masses. Moreover, another 
set of values of fi a and fis are currently being used at (3 — 1.90 to generate 
ensembles with a slightly lower D meson mass and a third value of the kaon 
mass, in order to properly interpolate the lattice data to the physical strange 
quark mass. 
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Fig. 3. (a): 2m 2 K — rrip S , and (b): mp, as a function of ntpg, for j3 = 1.95 (blue) 
and P = 1.90 (orange). The physical point is shown (black star). The kaon and D 
meson masses appear to be properly tuned at f3 = 1.95. The ensembles at (3 = 1.90, 
lis = 0.190 have a larger value of the strange quark mass, while the red point at 
(3 = 1.90, cilis = 0.197 appears to be well tuned. Data points have been scaled with 
the lattice spacing a = 0.08585(53) fm for /3 = 1.90, and a = 0.07820(59) fm for 
(3 = 1.95, obtained in this work and where the errors are only statistical. 



Ensemble 


P 


arrix 


arri£) 


A30.32 


1.90 


0.25150(29) 


0.9230(440) 


A40.32 




0.25666(23) 


0.9216(109) 


A40.24 




0.25884(43) 


0.9375(128) 


A40.20 




0.26130(135) 


0.8701(152) 


A50.32 




0.26225(38) 


0.9348(173) 


A60.24 




0.26695(52) 


0.9298(118) 


A80.24 




0.27706(61) 


0.9319(94) 


A100.24 




0.28807(34) 


0.9427(99) 


Al00.24s 




0.26502(90) 


0.9742(133) 


B25.32 


1.95 


0.21240(50) 


0.8395(109) 


B35.32 




0.21840(28) 


0.8286(85) 


B55.32 




0.22799(34) 


0.8532(62) 


B75.32 




0.23753(32) 


0.8361(127) 


B85.24 




0.24476(44) 


0.8650(76) 



Table 4 

For each ensemble, the values of the kaon mass and the D meson mass as determined 
in [66]. 

3.4 Discretisation Effects in Light-quark Observables 



In this section we explore discretisation effects in the analysed light-quark 
observables. To this aim we also make use of the determination of the chi- 
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are given in tables 1 and 9. 
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Fig. 5. The ratio "ip s //p S as a function of mpg/m^-, for the N{ = 2+1 + 1 ensembles 
at P = 1.90 and /3 = 1.95, compared to the N { = 2 data at /3 = 3.90, /3 = 4.05 and 
13 = 4.20 [17]. The physical point is shown (black star). 

rally extrapolated r value for our data samples, as discussed in the following 
section 3.5. 

In figures 4(a) and 4(b) we study the sensitivity of the charged pion mass 
and decay constant to possible discretisation effects, by comparing the Nf = 
2 + 1 + 1 data at (3 = 1.90 and = 1.95 and the results obtained in twisted 
mass simulations with two dynamical flavours [17]. The alignment of all data 



14 




2e-05 4e-05 6e-05 8e-05 1e-04 2e-05 4e-05 6e-05 
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Fig. 6. The Sommer scale r^/a as a function of (ani) 2 for (a) (3 = 1.90 and (b) 
(3 = 1.95. The lines represent a linear extrapolation in (a/i/) 2 to the chiral limit. 
The lightest point (open symbol) is not included in the fits and we have always used 
the largest available volume for a given value of the mass. 

points at different values of /3 is in itself an indication of small discretisation 
effects. The comparison and good agreement with the Nf = 2 data seems also 
to suggest no significant dependence upon the inclusion of dynamical strange 
and charm quarks for these light observables, at least at the present level of 
accuracy and provided that no cancellations occur due to lattice discretisation 
effects. However, only a more complete study at significantly different lattice 
spacings will allow to draw conclusions. 

In the same spirit, we show in figure 5 an analogous ratio plot where the nu- 
cleon mass data points are included. The alignment of all data and the good 
extrapolation to the physical point is again evident. We defer to future pub- 
lications the analysis of the baryon spectrum and the study of discretisation 
effects in strange- and charm-quark observables. 



3.5 The Sommer Scale r 

The Sommer scale r [67] is a purely gluonic quantity extracted from the static 
inter-quark potential. Since the knowledge of its physical value remains rather 
imprecise, we use the chirally extrapolated lattice data for r^/a only as an 
effective way to compare results from different values of the lattice spacing. 
In this work, the lattice scale is extracted by performing x?T inspired fits to 
the very precise data for af PS and am PS , and by using the physical values of 
m n and f n as inputs. 
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Figures 6(a) and 6(b) display the data for r^/a at both values of the lat- 
tice coupling = 1.90 and 1.95, and as a function of the bare lattice mass 
squared. The data are reasonably well described by a quadratic dependence, 
as also previously found for our Nf = 2 ensembles. For a more detailed dis- 
cussion of the possible functional forms and their theoretical interpretation 
see [37]. To extrapolate to the chiral limit, we have performed fits using the 
largest available volume at each value of the pseudoscalar mass. The chirally 
extrapolated values for our Nf — 2 + 1 + 1 ensembles are rj/a = 5.231(38) at 
f3 = 1.90 and r^/a = 5.710(41) at f3 = 1.95, where the lightest points of both 
ensembles have been excluded from the extrapolation, consistently with the 
fact that they do not satisfy our most stringent tuning condition to maximal 
twist. 

In order to meaningfully compare the dependence upon the light quark mass 
at the two different lattice couplings /3 = 1.90 and 1.95, we estimated the 
slope of the functional form tq/tq = 1 + c r (rg mps) 4 , where the explicit lattice 
spacing dependence has been removed. We observe a mild dependence on 
the light quark mass and similar slopes c r [(3 = 1.90] = —0.0379(37) and 
c r [(3 = 1.95] = —0.0234(69). It is also worth noticing that the dependence 
upon the light quark mass of the Nf — 2 + 1 + 1 data and that observed in 
the Nf = 2 case [37] are not significantly different. 

3. 6 Effects of Isospin Breaking 

A most delicate aspect of the twisted mass formulation is the breaking of 
the isospin symmetry. Clear evidence for this breaking has been found in the 
Nf = 2 simulations by ETMC when comparing the neutral with the charged 
pion masses. Indeed, while the discretisation effects in the charged pion were 
observed to be very small, significant 0(a 2 ) corrections appear when studying 
the scaling to the continuum limit of the neutral pion [17]. Notice, however, 
that similar effects have not been observed in other quantities that are in 
principle sensitive to isospin breaking but not trivially related to the neutral 
pion mass. These observations are supported by theoretical considerations 
detailed in [31,32]. 

In the Nf = 2 + 1 + 1 case, it turns out that the isospin breaking effect 
in the mass difference of charged and neutral pion masses is larger than for 
Nf = 2 at fixed physical situational as can be inferred from table 5. On the 
other hand, the same theoretical considerations as in [32] do apply to the 
case of Nf = 2 + 1 + 1 flavours, and it is expected that the same class of 
physical observables as for Nf = 2 will not be significantly affected by isospin 

3 Notice however that different gauge actions are used in the Nf = 2 and Nf = 
2 + 1 + 1 cases as described in section 2.1. 
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Ensemble 




r 0™PS 


y o 
r m PS 


c 


B35.32 
B55.32 


1.95 


0.7196(57) 
0.8861(67) 


0.388(40) 
0.679(40) 


-12.0(1.1) 
-10.6(1.8) 


B 6 N f = 2 
B 2 N f = 2 


3.90 


0.7113(66) 
0.9001(86) 


0.585(43) 
0.712(54) 


-4.6(1.5) 
-8.6(2.2) 



Table 5 



Measurements of the masses of the charged and the neutral pion. We compare 
runs at /3 = 1.95 and Nf = 2 runs [17] with comparable lattice spacing and similar 
charged pion masses in physical units. All masses are reported in units of the chir ally- 
extrapolated ro for the same ensemble, see table 9, and r^/a = 5.316(49) for Nf = 2. 
We also report on the approximate value of c, giving the slope of the a 2 dependence 
of the pion mass splitting. 

breaking corrections. Having said that, a careful measure of this effect for 
each observable or class of observables is anyway mandatory. The increase 
of the pion mass splitting with increasing the number of flavours in the sea 
is in line with the observation [35] of a stronger first order phase transition 
when moving from N{ — 2 to Nf — 2 + 1 + 1, as discussed in section 2.1. 
Indeed, the endpoint of the phase transition [44,45] corresponds to the critical 
value of the light twisted mass /i/ jC where the neutral pion mass vanishes. The 
mass difference can be described by ?"o 2 (( m Ps) 2 ~~ ( m ps) 2 ) = c ( a / r o) 2 > where 
the coefficient c is related to /^ iC [44, 45] and it is therefore a measure of the 
strength of the first order phase transition. Hence, a larger value of c means 
that simulations are to be performed at smaller values of the lattice spacing 
to reach, say, the physical point. Table 5 reports on the values of ffip S , mj s 
and c for some examples taken from the /3 — 1.95 ensemble and the Nf = 2 
ensemble with the closest values of the lattice spacing and physical charged 
pseudoscalar mass. As anticipated, the coefficient c increases in absolute value 
from N f = 2 to Nf = 2 + 1 + 1. 

We are currently performing simulations at a significantly different and lower 
lattice spacing than the present ensembles. They will allow to determine the 
slope c for Nf = 2 + 1 + 1 more accurately and to better quantify the conditions 
to approach the physical point. 

3. 7 Stout Smeared Runs 

In addition to our main simulation ensembles, we also performed runs with 
stout smeared gauge fields in the lattice fermionic action. The stout smearing 
as introduced in [59] was designed to have a smearing procedure which is 
analytic in the unsmeared link variables and hence well suited for HMC-type 
updating algorithms. In an earlier work with Nf = 2 quark flavours [60] we 
showed that using smeared gauge fields in the fermion operator is reducing 
the strength of the phase transition in twisted quark mass simulations and 
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Ensemble 


P 


^crit 


api 


ap a 


aps 


A tra j. 


ro/a 


A st 40.24 


1.90 


0.145512 


0.0040 


0.170 


0.185 


1500 


5.304(35) 


A st 60.24 




0.145511 


0.0060 






3100 


5.300(37) 


A st 80.24 




0.145510 


0.0080 






2000 


5.353(43) 



Table 6 

Parameters of the runs with stout smearing on L/a = 24, T/a = 48 lattices. The 
number of thermalised trajectories with length r = 1 is given by A^raj.- The label 
"st" in the ensemble name refers to the use of stout smearing, compared to the non 
stout-smeared ensemble in table 2. 



Ensemble 


ampg 


arriK 


arriD 




A st 40.24 
A st 60.24 
A st 80.24 


0.12600(93) 
0.14888(78) 
0.17156(69) 


0.2479(18) 
0.25338(67) 
0.26198(80) 


0.802(27) 
0.825(26) 
0.811(12) 


0.0175(68) 
0.0017(50) 
0.0138(48) 



Table 7 



The masses in lattice units for the ensembles with one level of stout smearing, 
therefore allows to reach smaller quark masses at a given lattice spacing. 

The definition of the stout smeared links can be found in [59] , and for the pa- 
rameter p connecting thin to fat gauge links we choose p = 0.15. In principle, 
such smearing can be iterated several times, with the price of rendering the 
fermion action delocalised over a larger lattice region. We made a conservative 
choice to maintain the action well localised and performed a single smearing 
step. As shown in [60], this kind of smearing does not substantially change 
the lattice spacing, and for the sake of comparison we thus kept the same 
value of f3 as in one of the non stout-smeared runs. On the other hand, the 
hopping parameter has to be tuned again, since the additive renormalisation 
of the quark mass is expected to be smaller. The parameters of our runs are 
given in Table 6. These runs have been done with the two-step polynomial Hy- 
brid Monte Carlo (TS-PHMC) update algorithm [68]. Results for the hadron 
masses are collected in Table 7, where the quoted errors include an estimate 
of the systematic error induced by variations of the fitting range. The method 
of estimating and combining statistical and systematic errors for the case of 
the kaon and D meson masses is described in [66]. 

As the values of mpcAc//-*/ m table 7 show, the hopping parameters are well 
tuned to maximal twist. The masses in the run with smallest light twisted 
mass api = 0.0040 (ensemble A st 40.24) satisfy r mpg = 0.668(10), r m^ = 
1.315(13) and r^mo = 4.25(29). This means that the pion is lighter than in 
the corresponding run without stout smearing (see table 8) and the kaon and 
D meson masses are closer to their physical value. The smaller pion mass 
should be interpreted as due to a quark mass renormalisation factor closer 
to one. For the same reason the tuned twisted masses in the heavy doublet 
ap a = 0.170, aps = 0.185 are smaller than in the runs without stout smear- 
ing. It is also interesting to compare the mass splitting of the charged and 
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neutral pion between runs with and without stout smearing. For the ensemble 
A st 60.24 we obtain a neutral pion mass r^mpg = 0.409(34) and a charged 
pion mass r^mpg = 0.7861(56), in units of the chirally extrapolated value 
rj/a = 5.280(25), providing an estimate of the slope c = —12.6(0.8). Notice 
that the mass dependence of r /a in table 6 is reduced as compared to the runs 
with no stout smearing, and a quadratic dependence on the bare quark mass 
has been used for the extrapolation to the chiral limit, consistently with the 
analysis of section 3.5. For the corresponding ensemble A60.24 without stout 
smearing, using data in tables 8 and 9, we obtain instead r^mpg = 0.560(37), 
r o m PS = 0.9036(71), and a slope c = —13.8(1.2), slightly but not significantly 
different from the stout-smeared case. 

The runs with stout-smeared gauge links show somewhat better character- 
istics than the ones without stout smearing, but the improvements are not 
dramatic, at least with one level of stout smearing. More iterations would 
further accelerate the approach to lighter masses and are expected to further 
reduce the charged to neutral pion splitting. However, it is a delicate matter to 
establish how physical observables other than the spectrum will be affected. 
Based on these considerations and given the present pool of data, the final 
results in this study are obtained with non stout-smeared simulations. 



4 Results: / PS , mps and Chiral Fits 

We concentrate in this section on the analysis of the simplest and phenomeno- 
logically relevant observables involving up and down valence quarks. These are 
the light charged pseudoscalar decay constant fp$ and the light charged pseu- 
doscalar mass m PS . 

The present simulations with dynamical strange and charm quarks, sitting at, 
or varying around, their nature given masses, should allow for a good measure 
of the impact of strange and charm dynamics on the low energy sector of 
QCD and the electroweak matrix elements. As a first step, one can determine 
the low energy constants of chiral perturbation theory (%PT). The values of 
afps and amps for our ensembles at (3 = 1.95 and (3 = 1.90 are summarised 
in table 8. In contrast to standard Wilson fermions, an exact lattice Ward 
identity for maximally twisted mass fermions allows for extracting the charged 
pseudoscalar decay constant fps from the relation 



without need to specify any renormalisation factor, since Zp = [3]. We 

have performed fits to NLO SU(2) continuum x?T at j3 — 1.95 and (3 = 1.90, 
separately and combined. Results are summarised in table 9. 
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Ensemble 


am 


ampg 


a/ps 


ro/a 


LI a 


A3n 39 


n nnsn 


n 1 9395(36)(1 4) 


n n645i (35)(3) 


5 91 7(3n) 




A4D 39 


n nn4n 


n 14149(97)(49) 


n nfi79i (i 


5 179(49) 


39 


A40.24 


0040 


n 1 4492(59) (34) 


n n6568(34)(7) 


5 178(44) 


24 


A40.20 


0040 


14871 (92)(11fi) 


n61 94(65) (23) 




20 


A 50 32 


0050 


15796(32)(28) 


07048(1 6") (4) 


5 081(45) 


32 


A60.24 


0060 


1 7275(45) (23) 


071 69(22)(2) 


5 909(58) 


24 


A80.24 


0080 


19875(41)(35) 


07623(21") (4) 


4 989(40) 


24 


A 1 nn 94 


n m nn 

\j . \J x \J\J 


n 22293(35)(38) 


n n7996(9n)(4) 


4 864(91 ) 


94 


A100 24s 


0.0100 


22125(58)(119) 


07843(26) (21) 


4 918(50) 


24 


B25.32 


0.0025 


0.10680(39)(27) 


0.05727(36)(8) 


5.728(35) 


32 


B35.32 


0.0035 


0.12602(30)(30) 


0.06074(18)(8) 


5.634(43) 


32 


B55.32 


0.0055 


0.15518(21)(33) 


0.06557(15)(5) 


5.662(33) 


32 


B75.32 


0.0075 


0.18020(27)(3) 


0.06895(17)(1) 


5.566(44) 


32 


B85.24 


0.0085 


0.19396(38)(54) 


0.06999(20)(5) 


5.493(41) 


24 



Table 8 



Lattice measurements of the charged pseudoscalar mass amps, the charged pseu- 
doscalar decay constant a/ps and the Sommer scale in lattice units ro/a for our two 
ensembles at f3 = 1.90 (A set) and (3 = 1.95 (B set). The value of the light twisted 
mass afii and the spatial length L/a are also shown. Quoted errors are given as 
(statistical) (systematic) , with the estimate of the systematic error coming from the 
uncertainty related to the fitting range. 

We thus simultaneously fit our data for the pseudoscalar mass and decay 
constant to the following formulae, where the contributions F, D and T 
parametrising finite size corrections, discretisation effects and NNLO x?T 
effects, respectively, will be discussed below: 

™ls(L) = Xm (l + ^3 + D m 2 s a 2 + e T mls ) F m%s 

/ps(L) = f (l - 2 ^4 + D fps a 2 + e TfJ) F fps , (11) 

with the pseudoscalar mass squared at tree level denned as Xn = 2 B fii and 
the chiral expansion parameter by £ = x^l (4tt/o) 2 - The low energy constants 
Is and U receive renormalization corrections according to k — k + In [A 2 /x^], 
with A the reference scale. During the fitting procedure, where all quantities 
are denned in lattice units, we set the reference scale to a single lattice spac- 
ing to let its constant logarithmic contribution vanish. Once the scale of the 
simulation has been set, the low energy constants are rescaled to the scale of 
the physical pion mass to recover the physical values 1$ and I4. 

Systematic errors can arise from several sources: finite volume effects, neglect- 
ing of higher orders in %PT and finite lattice spacing effects. These different 
corrections are accounted for explicitly in eq. (11). Finite volume corrections 
are described by the rescaling factors denoted by F m 2 s and -F/ PS , computed 
in the continuum theory. Notice that the discretisation effects present in the 
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neutral pion mass, see section 3.6, generate peculiar finite volume corrections 
which have been recently analysed in ref. [69]. We shall comment on them 
later. We investigated the effectiveness of one loop continuum finite vol- 
ume corrections, as first computed in [70], which do not introduce any addi- 
tional low energy constants. However, the resummed expressions derived by 
Colangelo, Diirr and Haefeli (CDH) in [71] describe the finite volume effects in 
our simulations better, be it at the expense of the introduction of two new free 
parameters, and are thus adopted for this analysis. To 0(£ 2 ), these corrections 
read 



F 2 

m PS 



oo 

l _Y_P*_U I ($ + e2 I (4) 
n=l 1 A n 



r(4) 
7 



with geometric contributions defined as 



(12) 



4 2 ) = -2i^(A n ) 

4 4) = - f * + &h + ~ l 2 ~ 5/ 3 -4/4) tfi(A„) + 

238 61 16 , 64 \ K 2 (K 
n — 7T h — 7rh 



9 6 3 3 
if = -AK 1 (X n ) 

4 4) =(|-|- + 4/ 1 + ^ 2 -6/ 4 ) K 1 (X n ) + 

307 391 16 64 \ K 2 (X n ) 

+ — — 7T — Li Trh) — : • (13) 



9 24 3 3 V A„ 

The ii'i are the modified Bessel functions and the low energy constants l\ and 
I2 again receive renormalisation corrections. Equations (12) and (13) use the 
shorthand notation A n = i/nmpgl. The p n in eq. (12) are a set of multiplic- 
ities, counting the number of ways n 2 can be distributed over three spatial 
directions^ Because the finite volume corrections in the case of the volumes 
used in the chiral fits are fairly small to begin with and subsequent terms 
quickly decrease, the sums over n can be truncated rather aggressively with- 
out real loss of precision. It is therefore unnecessary, in practise, to go beyond 
the lowest contributions. The parameters l\ and Z 2 , which are in fact low en- 
ergy constants appearing at NLO in xPT, cannot be determined well from the 
small finite volume corrections alone. Priors are therefore introduced as addi- 
tional contributions to the % 2 , weighting the deviation of the parameters from 

These values are straightforwardly precomputed to any order, but are also given 
in, e.g. [71]. 
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their phenomenological values by the uncertainties in the latter. The values 
used as priors are -0.4(6) for l\ and 4.3(1) for l 2 [71], as reported in table 9. 
We used the largest available volumes for each ensemble, in the ;^PT fits. For 
those points, the difference between the finite volume and the infinite volume 
values estimated via CDH formulae for /pg and mp s are within 1%, except for 
the runs B85.24 and A60.24 (see table 2 and table 8), where they are about 
1.5% for both quantities. 

Because of the automatic 0(a) improvement of the twisted mass action at 
maximal twist, the leading order discretisation artefacts in the chiral formulae 
of (11) are at least of 0(a 2 ), and 0(a 2 n) for mp s . The mass and decay constant 
of the charged pion have been studied up to NLO [44, 45, 50] in the context 
of twisted mass chiral perturbation theory (tm^PT). The regime of quark 
masses and lattice spacings at which we have performed the simulations is 
such that [Mi > aAq CD . In the associated power counting, at maximal twist, 
the NLO tm^PT expressions for the charged pion mass and decay constant 
preserve their continuum form. The inclusion of the terms proportional to 
D m 2 f FS , parametrising the lattice artifacts in eq. (11), represents an effective 
way of including sub- leading discretisation effects appearing at NNLO. The 
finite lattice spacing artefacts can of course not be determined using only data 
from a single lattice spacing. In addition, including these terms when analysing 
data with an insufficient range in a, may lead to mixing of these degrees of 
freedom with continuum parameters and thereby destabilise the fits. Hence, 
these terms were neglected for the separate fits, but included to arrive at a 
qualitative estimate of these systematic effects in a combined fit to the data 
at both lattice spacings. 

Finite size effects on our data at finite lattice spacing can be analysed in the 
context of twisted mass chiral perturbation theory as recently proposed in 
ref. [6 9], p] However, our present limited set of data with only a small number 
of different volumes all of them at a single value of the lattice spacing, is not 
sufficient to apply such an analysis. We plan, however, to perform dedicated 
runs on different volumes to confront our data to the finite size effect formulae 
of ref. [69] and to estimate in particular the size of the pion mass splitting in 
this alternative way. 

Finally, results from continuum ;\PT at NNLO can be included to examine 
the effect of the truncation at NLO. They are given by 



5 Notice that, in principle, after performing the continuum limit at fixed physical 
volume, finite size effects can be analysed by means of continuum xPT. 
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T mls = — (49 + 28 Z! + 32/ 2 - 9/ 3 ) + 4 k m 

T /ps = ~ (23 + 14 h + 16 / 2 + 6 h - 6 Z 4 ) + 4 (14) 

Two new parameters fc m and fc/ enter these corrections. Again, a limited range 
of input pion masses may lead to poorly constrained values of these newly 
introduced parameters, some degree of mixing among different orders and fit 
instabilities. To retain predictive power and stability, additional priors are 
given for k m and kf, both priors set to 0(1), analogously to what is done for 
li and I2 in the CDH finite volume corrections. 

To set the scale at each lattice spacing, we determine o/Xphys; 

the value of a\ii 

at which the ratio \Jm^ s {L = oo)// P g(L = 00) assumes its physical value. We 
can then use the value of / P g, or equivalently mps, to calculate the lattice 
spacing a in fm from the corresponding physical value. We also perform a 
chiral fit combining the two different lattice spacings. With only two different 
values of j3, that are in fact fairly close to each other, a proper continuum 
limit analysis cannot be performed. Instead, we treat this combined fit as a 
check on the presence of lattice artefacts and the overall consistency of the 
data. Without a scaling variable, such as the Sommer scale tq, the data from 
different lattice spacings cannot be directly combined. Rather, the ratios of 
lattice spacings and light quark mass renormalisation constants (Z^ = 1/Z P ), 
as well as the renormalised B parameter are left free in the fit. 

In order to estimate the statistical errors affecting our fitted parameters, we 
generate at each of the \i\ values 1000 bootstrap samples for m PS and / PS 
extracted from the bare correlators, organised by blocks. For each sample, and 
combining all masses, we fit mp s and / P g simultaneously as a function of \x\. 
The parameter set from each of these fits is then a separate bootstrap sample 
for the purposes of determining the error on our fit results. By resampling /pg 
and m PS on a per-configuration basis, correlations between these quantities 
are taken into account. 

Our final results for the separate and combined fits are summarised in table 9. 
The ;\TT fit ansatze provide a satisfactory description of the lattice data, with 
a xVd.o.f = 5.68/3 ~ 1.9 at /3 = 1.95, x 2 /d.o.f = 4.31/5 ~ 0.9 at (3 = 1.90, 
and 16.9/11 ~ 1.5 for the combined fit. We also predict the scalar radius of 
the pion at next to leading order 

< r2 >- L0 = ('< " n) • (15) 

The numerical values in table 9 for the combined fit show a very good agree- 
ment with the results from the separate fits, and with errors at the percent 
level throughout. The fits for / PS and m PS at /3 — 1.95 are displayed in fig- 
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/3 = 1.90 


/3 = 1.95 


combined 


priors 


^3 


3.435(61) 


3.698(73) 


3.537(47) 


- 


k 


4.773(21) 


4.673(25) 


4.735(17) 


- 


h 


-0.296(104) 


-0.430(93) 


-0.309(139) 


-0.4(6) 


h 


4.260(12) 


4.329(15) 


4.325(10) 


4.3(1) 


fo [MeV] 


120.956(70) 


121.144(83) 


121.031(54) 




fn/fo 


1.0781(18) 


1.0764(18) 


1.0774(17) 




on / 2 


l. 029(16) 


l. 032(21) 


l. 030(13) 




/ 2\NT,0 rr 2i 

(r)s u [fm J 


0.7462(43) 


0.7237(51) 


0.7375(34) 




X 1 ( o -i r\r\\ 

rg/a(p = 1.90) 


5.231(38) 




5.231(37) 




r$/a(l3 = 1.95) 




5.710(41) 


5.710(42) 




rJ(/3 = 1.90) [fm] 


0.4491(43) 




0.4505(40) 




rJ(/3 = 1.95) [fm] 




0.4465(48) 


0.4439(39) 




a(/3 = 1.90) [fm] 


0.08585(53) 




0.08612(42) 




a(/3 = 1.95) [fm] 




0.07820(59) 


0.07775(39) 





Table 9 

Results of the fits to SU(2) xPT for the ensembles at f3 = 1.95 and (3 = 1.90, 
separate and combined. The largest available volumes are used for each ensemble. 
Predicted quantities are: the low energy constants 1%^ (while l\ : 2 are introduced 
with priors), the charged pseudoscalar decay constant in the chiral limit fo, the 
mass ratio 2£?o/U//mp S at the physical point and the pion scalar radius (r 2 )^ L °. 

The scale is set by fixing the ratio /p^ _ °°' /mp^ -00 ' = fir/m^ = 130.4(2)/135.0 to its 
physical value [41]. The chirally extrapolated Sommer parameter is determined 
separately and not included in the chiral fits. For a comparison with the Nf = 2 
ETMC results, see [17]. 

ures 1(a) and (b), while in figures 7(a) and (b) we show the analogous fits at 
(5 = 1.90. Figures 8(a) and (b) show the results for the fit combining the two 
/3 values. 

The data presented here do not allow yet for a complete account of the system- 
atic effects, but we extract estimates of their magnitude by extending the fits 
with additional terms as written down in eq. (11). Checks were done for ;yPT 
NNLO terms and 0(a 2 ) corrections separately. Including NNLO corrections 
does not lower the total x 2 °f the fit, while we do observe a shift of several 
standard deviations for the lower order parameters already present in the NLO 
fit. Using these shifted values to obtain the implied NLO approximation pro- 
duces fits with much larger values of \ 2 ■ We conclude that the current data 
lack the precision and range in quark masses to constrain NNLO effects, the 
added degrees of freedom mix with NLO effects and destabilise the fit instead. 
In practise, we conclude that the systematic error from the truncation of xPT 
is unobservable at the current level of precision. Inclusion of 0(a 2 ) corrections 
leads to similar observations, as the difference between the lattice spacings 
and the statistical accuracy of the data is too small to result in a stable fit. 
The fit mixes Df PS and D m 2^ on the one hand and f , B and the rescaling in 
the lattice spacing and the quark mass on the other. 
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(a) (b) 

Fig. 7. (a) The charged pseudoscalar mass ratio rrip S /2Bofii and (b) the pseu- 
doscalar decay constant fps as a function of 2Bq[ii, for the ensemble at /3 = 1.90, 
fitted to SU(2) chiral perturbation theory, eq. (11). The scale is set by a// p hy S > 
the value of apii at which the ratio /Ji^ - ^ / m ps _<X ^ assumes its physical value [41] 
fir/ m ir = 130.4(2)/135.0 (black star). The light twisted masses used in the fit range 
from a\i[ = 0.004 to 0.010. The lightest point (open symbol) lies outside our most 
conservative tuning criterion to maximal twist, and is not included in the fit. 

The chirally extrapolated Sommer scale Tq has been determined separately, 
using a fit of r^/a with quadratic dependence on the bare light quark mass, 
as shown in figures 6(a) and 6(b), and using the lattice spacing determined 
by the chiral fits. As also reported in table 9, the obtained values are = 
0.4491(43) fm at p = 1.90 and = 0.4465(48) fm at p = 1.95, where only 
statistical errors are quoted. For consistency, we also verified that a combined 
chiral fit with the inclusion of Tq/cl, as data points and additional fit parameter, 
gives results anyway in agreement with the strategy adopted here. 

For our final estimates of the low energy constants ^4 and the chiral value of 
the pseudoscalar decay constant f we use the predictions from the (3 = 1.95 
ensemble based on two important observations. First, the strange quark mass 
in this ensemble is better tuned to the physical value. Secondly a reduced 
isospin breaking is observed at this finer lattice spacing. The results for the 
P = 1.90 ensemble and the combined fits serve instead as an estimation of 
systematic uncertainties. As a result of the current Nf = 2 + 1 + 1 simulations 
we thus quote 

l 3 = 3.70(7) (26) T A = 4.67(3) (10), (16) 

and fo = 121.14(8)(19) MeV, where the first error comes from the chiral fit 
at P = 1.95, while the second quoted error conservatively accommodates the 
central values from the P = 1.90 and combined fits as a systematic uncertainty. 
The predictions for l 3 and Z4 are in good agreement and with our two-flavour 
predictions [17] and with other recent lattice determinations [2,72]. 
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Fig. 8. (a) The charged pseudoscalar mass ratio (mps/2i?ow) 2 an d 0°) the pseu- 
doscalar decay constant /ps as a function of 2£>o/^, for the combined ensembles at 
jS = 1.90 and j3 = 1.95, and fitted to eq. (11). The scale is set as in figure 7 (black 
star). The light twisted masses used in the fit range from afii = 0.0035 to 0.010. 
The lightest point at f3 = 1.90 (open orange symbol) and at f3 = 1.95 (open blue 
symbol) lie outside our most conservative tuning criterion to maximal twist, and 
are not included in the fit. 

5 Conclusions and Outlook 

In this paper we have presented the first results of lattice QCD simulations 
with mass-degenerate up, down and mass-split strange and charm dynami- 
cal quarks using Wilson twisted mass fermions at maximal twist. This study 
constitutes a first step in our effort to describe low energy strong dynamics 
and electroweak matrix elements by fully taking into account the effects of a 
strange and a charm quark. 

We have considered ensembles at slightly different lattice spacings simulated 
with Iwasaki gauge action at (3 = 1.95 with a pa 0.078 fm and /3 = 1.90 with 
a pa 0.086 fm. The charged pseudoscalar masses range from 270 to 510 MeV 
and we performed fits to SU(2) chiral perturbation theory with all data at 
a value of m-psL > 4. This analysis provides a prediction for the low energy 
constants I3 = 3.70(7)(26) and I4 = 4.67(3)(10), for the charged pseudoscalar 
decay constant in the chiral limit / = 121. 14(8) (19) MeV and for the scalar 
radius at next-to-leading order (r 2 )^ LO = 0.724(5) (23) fm 2 . A companion 
paper [66] describes the less straightforward determination of the kaon and 
D-meson masses for the same ensembles. 

We have compared our results in the light meson sector with those obtained 
for iVf = 2 flavours of maximally twisted mass fermions, ref. [17]. There, 
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an extrapolation to the continuum limit, a study of finite size effects and 
checks against higher order xPT have been performed, leading to a controlled 
determination of systematic errors. The comparison we have carried through 
does not show any significant difference between N { = 2 and Nf — 2 + 1 + 1 
flavours, at least at the present level of accuracy. These results would suggest 
that effects of the strange and charm quarks are suppressed for these light 
observables, as it should be expected. The same comparison has also been 
used for a first investigation of lattice discretisation errors. As figures 4(a) 
and 4(b) show, the Nf = 2 + f + 1 data are completely consistent with the 
corresponding ones obtained for Nf = 2, where the discretisation effects have 
turned out to be very small. Thus, it can be expected that also for the case 
of Nf = 2 + 1 + 1 flavours the lattice spacing effects will be small, at least for 
the light meson sector considered here. Notice however that, at the present 
level of accuracy, there is still the possibility that cancellations occur between 
physical contributions due to dynamical strange and charm quarks and lattice 
discretisation effects. A more accurate study at a significantly lower lattice 
spacing will allow to draw conclusions. 



One aspect of the twisted mass formulation is the breaking of isospin sym- 
metry. Its effect is likely to be most pronounced in the lightest sector, where 
lattice discretisation effects at 0(a 2 ), affecting the neutral pseudoscalar mass 
only, generate a mass splitting between the charged and the neutral pseu- 
doscalar mesons. While this mass splitting for Nf = 2 + 1 + 1 flavours has 
been found here to be larger than in the Nf = 2 simulations at fixed physical 
situation, we do not find further effects in other quantities computed so far. 
This observation is supported by theoretical arguments [31,32] and consistent 
with our experience in the Nf = 2 flavour case. 



We consider the present results to be encouraging to proceed with the Nf = 
2 + 1 + 1 flavour research programme of ETMC. In particular, we want to 
perform the non-perturbative renormalisation with dedicated runs for Nf = 4 
mass-degenerate flavours, an activity which we have started already. Further- 
more, we want to compute the quark mass dependence of many physical quan- 
tities towards the physical point where the pion assumes its experimentally 
measured value. We are currently performing simulations at a significantly dif- 
ferent and lower lattice spacing than the present ensembles. Both strategies, 
smaller quark masses and smaller lattice spacings, will allow us to estimate 
systematic effects on a quantitative level and to obtain in this way accurate 
physical results in our Nf = 2 + 1 + 1 flavour simulations with statistical and 
systematical errors fully under control. 
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